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circular  orbit 
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Hamilton's  canonical  equations.  An  algorithm  Is  then 
developed  to  find  closed  orbits  for  a  specified  period. 
Periodic  orbits  in  Phobos*  orbital  plane  have  been  found  to 
be  unstable.  A  control  system  is  developed  using  scalar 
modal  control  to  stabilize  the  orbit  by  shifting  the 
unstable  Poincar*  exponent  to  the  left-half  plane. 
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Introduction 


I . 


There  are  several  moons  in  the  solar  system  which  may 
be  modeled  as  triaxial  ellipsoids.  This  investigation 
studies  periodic  orbits  about  one  such  ellipsoid,  Phobos. 

Hamilton's  canonical  equations  are  used  to  derive  the 
equations  of  motion  for  a  satellite  in  orbit  around  Phobos. 
A  method  will  be  developed  to  locate  periodic  orbits  in 
Phobos'  orbital  plane,  and  determine  their  stability. 
Periodic  orbits  in  Phobos'  orbital  plane  have  been  shown  to 
be  unstable  (6:  30),  so  modal  control  theory  will  be  used  to 
develop  a  control  system  which  will  stabilize  these  orbits. 
Once  the  control  has  been  implemented,  it  will  be  tested  to 
determine  how  far  the  satellite  may  deviate  from  its  nominal 
trajectory  before  the  control  becomes  unstable.  These 
results  are  presented  in  graphical  format. 


Phobos  may  be  modeled  as  a  triaxial  ellipsoid.  (See 


Figure  1.)  The  coordinate  system  shown  in  the  figure 
coincides  with  Phobos'  center  of  mass.  Phobos  rotates 
about  it's  smallest  axis  (z  -  out  of  page),  with  its  longest 
axis,  y,  pointed  towards  Mars.  Therefore,  Phobos'  rotation 
rate  about  its  own  axis,  O,  equals  its  rotation  rate  about 
Mars,  w  (6:2).  For  this  investigation,  Phobos'  orbit  about 
Mars  is  taken  to  be  circular  with  radial  distance  jDj.  The 
satellite's  orbit  about  Phobos  is  denoted  by  R.  Appendix  A 
contains  the  needed  parameters  for  Mars  and  Phobos. 


Figure  1  :  Coordinate  System  for  Mars  and  Phobos 


Hamilton's  canonical  equations  vill  be  used  to  derive 


the  equations  o£  motion.  The  Hamiltonian  is  defined  as 

H  =  J  p.qt  -  L  (2-1) 

where  q^  and  are  the  generalized  coordinates  and  momenta, 
respectively.  L  is  the  Lagrangian  and  is  defined  as 


L  =  T  -  V  (2-2) 

where  T  is  the  satellite's  kinetic  energy  and  V  is  its 
potential  energy.  L,  and  therefore  T  and  V,  must  be  derived 
in  order  to  assemble  the  Hamiltonian. 

Kinetic  Energy  The  kinetic  energy  of  a  satellite  with 
mass,  is  defined  to  be 

1  * 

T  =  ^  M.«t  r.ot  (2-3) 

where  Tsai  is  the  velocity  of  the  satellite  measured  with 
respect  to  an  inertial  frame.  Using  the  center  of  Mars  as 
the  Inertial  reference  results  in 

z 

Vm*  t  =  d*d  (2-4) 

where  d  is  the  time  rate  of  change  of  d  with  respect  to  the 
inertial  frame.  Using  Figure  1,  d  may  be  written  as 


3 


u 


■ 


tz 


d  *  R  -  6  (2-5) 

where  D  is  the  time  derivative  of  D  with  respect  to  the 
inertial  frame  and  may  be  expressed  as 

6  =  Dr«l  ♦  w  x  D  (2-6) 

The  first  term  in  Equation  (2-6)  is  the  time  derivative  of  D 
with  respect  to  the  rotating  trame.  Since  Phobos  is  in  a 
circular  orbit  about  Mars,  this  term  is  equal  to  zero. 
Therefore, 


D  =  w  X  D  =  wz  X  Dy  *  -toDx 


Since  w  and  Q  are  equal 


D  =  — ODx 


(2-7) 


( 2—8 ) 


Likewise,  R,  the  time  derivative  of  R  with  respect  to  the 
inertial  frame,  may  be  written 


R  =  Rr*l  +  n  X  R 


If  R  is  defined  to  be 


R  =  Xx  +  Yy  +  Z* 


(2-9) 


(2-10) 


1 


4 


then 


■ 


r~ 


1 


R  =  (  Xx  ♦  Yy  +  Zz )  *  (flXy  -  OYx)  (2-11) 

Combining  Equations  (2-5),  (2—8),  and  (2—11)  results  In 

d  =  (X  -  OY  OD)x  +  (Y  ♦  Ctt)y  +  Zz  (2-12) 

Nov,  the  expression  for  kinetic  energy 

T  =  |  Moot  d*d  (2-13) 

may  be  written  using  Equation  (2-12)  as 

T  =  ^Maat  £(X  -  OY  +  OD)*  +  (Y  +  OX)*  +  Z*J  (2-14) 

Potential  Energy  The  potential  energy  of  the  satellite 
is  the  sum  of  the  potential  energies  due  to  Mars  and  to 
Phobos 


1 


Vaat  =  Vmors  ♦  VpKob  (  2—1 5 ) 

Potential  Energy  Due  to  Mars  The  potential  energy 
of  the  satellite  due  to  the  gravitational  attraction  of  Mars 
is 

„  GMn»«r»M»al  ,  _  ,  ,  . 

Vma rm  =  —  3 (  2—16  ) 

a 


5 


The 


where  G  is  the  Universal  Gravitational  Constant, 
magnitude,  d,  may  be  written 


d 


|R-Df 


2R*  D  + 


(2-17) 


the  binomial  expansion  of  which  may  be  shown  to  be 


d’1 


1 


3  f2DY  -  R2 

8l  D2 


(2-18) 


Since  D  »  R,  the  higher  order  terms  of  the  expansion  are 
negligible  and  may  be  Ignored.  Substituting  Equation  (2—18) 
in  Equation  (2—16)  provides  an  expression  for  the  potential 
due  to  Mars 


Vmor» 


(2-19) 


Potential  Energy  Due  to  Phobos  The  gravitational 
potential  for  an  arbitrary  body  is  derived  by  Meirovitch 
(5:  430-436)  and  will  be  used  to  develop  an  expression  for 
the  satellites  potential  due  to  Phobos.  Figure  2  shows  the 
coordinate  system  that  will  be  used  for  this  development. 
The  integral  form  of  gravitational  potential  energy  is 


Vphob  =  —  GMaot 


I 


dMpbob 


(2-20) 


6 


where 


2R*p  + 


(2-21) 


Doing  another  binomial  expansion  results  in 


-» 

r 


R 


-i 


+ 


R IE 

R8 


liBifiL*.  _  of  e3  I 

2RS  R5  J 


(2-2  2) 


7 


The  terms  of  Order 


and  above  in  Equation 


(2-22)  will 


be  dropped  so  that  the  equations  of  motion  will  be  a  system 
of  linear  differential  equations.  Using  the  following 
definitions  for  p  and  R 


—  A 

p  =  ax  + 

by  + 

A 

CZ 

(2 

50 

n 

X 

x> 

Yy  + 

z£ 

(2 

and  after 

some 

algebraic 

manipulation. 

Equations  ( 

(2—20)  may 

be  combined  to  give 

Vp**  .  GM-“  | 

d Mphob  — 

GM*at 

R* 

|  £aX  +  bY 

+  czj  dMphob 

iMi«t  f  r  3_ 

2RS  J  L  R2 


GM.at  |  |  J  (ax  +  by  +  cZ)2  -  (a2  +  b2  +  c2)  I  dMphob 


( 2-24) 


Evaluating  the  first  integral  in  Equation  (2—24)  gives 


GM«<u  f 

R  J 


dMphob  =  — 


GM»at  Mphob 


( 2-25) 


The  second  integral  in  Equation  (2-24)  goes  to  zero 
when  the  origin  of  the  axis  system  is  placed  at  the  center 
of  mass. 

The  third  integral  in  Equation  (2—24)  may  be  expanded 


and 


8 


j  (ax  +  bY  +  cZ)  -  (a  ♦  b  +  c  )  J  dMphob 


+ 

f  — -ll 

f  C*dMphob  +  —  2— 

I  ab  dMphob 

U*  J 

J  R2  J 

1 

+  f  ac  dMphob  +  $22  f  be  dMphob  (2-26) 

R2  J  R2  J 


By  using  definitions  for  aass  moments  of  inertia  (Ixx,  Iyy, 
Izz)  and  mass  products  of  inertia  (Ixy,  I xz,  Iyz),  the  terms 
of  Equation  (2-26)  may  be  simplified 


|  a2  dMphob  *  j  |  £(a2  +  b2)  +  ( 


a2  +  c2) 


-  (b 


2  ♦  c2)] 


dMphob 


(I**  ♦  Iyy  —  Ixx) 


J 


b  dMphob 


*  |  |  [(a*  +  b8)  +  <b*  +  °X)  ’  <a*  +  Ct)]  dMph°b 


j  (I**  +  I xx  —  Iyy) 


J*  C*  dMphob  *  J  J  |^( 


dMphob  *  7?  |  |  (a2  +  c2)  +  ( b2  +  c2)  -  (a2  +  b2)  |  dMphob 


1 

2 


(Iyy  + 


I  XX 


lx*  ) 


l 

I 

f 


ab  dMphob 


ac  dMphob 


be  dMphob 


Ixy 


Ixx 


Iy* 


(2-27) 


If  the  axes  are  arranged  such  that  they  align  with  the 
principle  axes  of  Phobos,  then  the  products  of  inertia  are 
zero.  Equation  (2—26)  may  now  be  rewritten 
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Vphob 


GMaat Mphob 
R 


GM»a t  Mphob 

4R* 


(Iyy  4  In  -  Ixx) 


4 


) 


(  Ixx 


4 


Itl 


Iyy) 


„.2 

■ 

p--l 

j ( Ixx  4  Iyy  - 

In) 

l  R* 

(2-28) 


which  further  simplifies  to 


Vphob 


GMaat  Mphob  GMaat  .  _  _ 

- - -  +  - —  (Ixx  4  Iyy  4 

R  4R* 


Ixx) 


--M~—  T  XZ(Iyy  4  In  -  Ixx)  4 

4R9  L 


Y*( Ixx 


+  I  **  -  Iyy)  +  Z  ( I  ** 


j 


(2-29) 


Hamiltonian  Mechanics  Since  expressions  for  both  the 
kinetic  and  potential  energies  have  been  developed,  the 
Lagrangian  may  be  formed.  For  convenience,  the  mass  of  the 
satellite  will  be  divided  out,  and  the  equations  developed 
will  be  on  a  per  unit  mass  basis.  Substituting  Equations 
(2—29),  (2—19),  and  (2-14)  into  Equation  (2—2)  gives  a 
Lagrangian 


L 


H 


(X  +  Q{ D  -  Y> )  +  (Y  ♦  OX) 


*  *  i*  ] 


GMmara 


[,  ,  DY  R*  3  f2DY  -  R*]  1 

D*  "  2D*  '  Sl  D*  J  J 


GMphob  G  , T  .  r  .  t  \ 

+  - = - - -  ( I xx  +  Iyy  +  I**) 

R  4R* 


^  [  ** 
4RS  «■ 


( I yy  +  I xs  —  I xx ) 


+  Y2(  Ixx  ♦  lx*  -  Iyy)  +  2Z(Ixx  +  Iyy  —  I**)  1  (2—30) 


Before  the  Hamiltonian 


H  *  PxX  +  PyY  +  PxZ  —  L 


(2-31) 


may  be  formed,  expressions  for  the  generalized  velocities 
are  needed.  This  Is  accomplished  by  using  Lagrange's 
equations 


at, 

Px  =  —  =  X  +  0(D  -  Y) 

ax 


( 2— 31a ) 


Py  =  —■  =  Y  +  OX 

ax 


( 2— 31b ) 


P*  =  =  z 

az 


( 2— 31c ) 


12 

f 


which  allow  the  Lagranglan  to  be  siapllfied  to 


L  *  \ 

£  Px*  +  1 

?y*  P.z  ] 

/  «\2  “I 

GMttxxtx  1 

1  +  DY 

R*  _ 

3  2DY  -  R  | 

*  D 

1  4-  -  — 

L  *>* 

2D* 

'  H  o*  J  J 

GMphob 

R 

r* 

-  -  (  I XX  +  Iyy  ♦  I**) 

4R* 

♦  1 
4R*  1 

X*(Iyy  +  1**  —  Ixx)  +  ¥*(Ixx  +  I**  - 

Iyy) 

+  Z*(Ixx  +  Iyy  —  Ixx)  j 

(2-33) 

This  allows  Equation  (2-31)  to  be  rewritten  as 

H  =  Pxjpx  + 

0(Y  -  D)j  +  Py[W  -  Oxj  +  Px*  -  L 

(2-34) 

The  Hamiltonian  may  now  be  presented  In  Its  final  fora 


where 

motion 


Px*  +  Py*  +  P.*  |  +  PxO(Y  -  D)  -  PyOX 


GHm«ur* 

,  x  DY 

R* 

3 

l'' 

2DY  -  R 

^2  - 

D 

l1*? 

1  n 

a 

1  CM 

1 

8 

.  D*  . 

GMphob 


[  X*(I 


2Ixx)  ♦ 


Y*( I  -  2Iyy)  +  Z2(I  -  21.*) 


] 


(2-35) 


I  =  I  xx  +  Iyy  ♦  I**.  Finally,  the  equations  of 

may  be  formed 


X 


OH 

«Px 


=*  Px  +  0<Y  -  D) 


(  2— 36a  ) 


a 

Y 


0H 

apx 


Py  -  nx 


(  2— 36b ) 


dw 


z 


*  p. 


(2-36c) 


Px  *  —  ^  =  QPy 


_  GMmara  j  X  _  3XY 
D  l  D*  D* 

^[<5X*  - 

9  4R  [ 


3GIX 

4R 


3R*X 

2D4 


2R* )  ( I 


GMphobX 


-  2Ixx) 


+  5Y*( I  -  lyy)  +  5Z*(I  -  I**) 


Py  ~  ~ 


2”  =  _OPx  _  GMmor. 
<*Y  *  D 


GMphobY  3GIY 


1  4Y  3 ( YR*—  DR*-  2DY*) 
D  D* 


2D 


4R* 


-  1Gy[  5x.(1  _ 

4R7  L 


2Ixx) 


+  ( 5Y*—  2R* ) ( I  -  lyy)  +  5Z*(I  -  la*) 


P* 


OH  _  GMwiar*  j  Z  _  3YZ  _  3R*Z 
*Z  ^  D  [  D*  D*  2D4  . 


GHphobZ 


mi.  2gz[5X*(I  _ 

4RS  4R? 


2 1  xx ) 


♦  5Y*(  I  -  lyy)  +  (  5Z*-  2R*)(I  -  Izz) 


These  are  the  equations  of  motion  for  this  system, 
also  listed  In  Appendix  B  for  convenience. 


( 2— 36d ) 


( 2— 36e ) 


( 2— 36f  ) 


They  are 
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III.  Periodic  Orbits 


The  equations  o£  motion  may  now  be  integrated  with 
selected  initial  conditions  in  an  attempt  to  locate  periodic 
orbits.  If  the  initial  conditions  do  not  produce  a  closed 
orbit,  the  equations  are  reintegrated  with  new  conditions 
until  a  closed  orbit  is  found.  To  prevent  arbitrary 
selection  of  initial  conditions,  an  algorithm  is  required 
which  will  systematically  alter  the  initial  conditions  until 
the  desired  orbit  is  found.  But,  before  the  algorithm  is 
developed,  it  would  be  beneficial  to  know  in  what  way  the 
orbit  trajectory  is  altered  when  the  initial  conditions  are 
changed.  There  exists  a  method  (8:  112-114)  to  obtain  this 
information  which  will  be  described  here  for  completeness. 

State  Transition  Matrix 

The  equations  of  motion  may  be  represented  as 

X  =  f(X)  (3-1) 

where  X  and  f(X)  are  vectors  whose  elements  are 


fx 

Eqn 

2-36a' 

Y 

Eqn 

2— 36b 

Z 

Px 

f(X)  =  ■ 

Eqn 

Eqn 

2— 36c 
2— 36d 

Py 

Eqn 

2— 36e 

Px 

Eqn 

2-  36  f 
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If  a  set  of  Initial  conditions  result  in  a  trajectory. 


Xo( t ) ,  then  a  nearby  trajectory  may  be  represented  by 

X(t)  *  Xo( t )  ♦  6X(t)  (3-3) 

where  <5X  is  a  small  deviation  from  the  initial  trajectory. 
Substituting  this  trajectory  into  the  equations  of  motion 
results  in 


X(t)  *  Xo( t )  ♦  <5X(t)  *  f(X  +  6X)  (3-4) 

Expanding  Equation  (3—4)  in  a  Taylor  series  about  the 
initial  trajectory  (ie,  about  6X  =  0)  produces 

Xo  +  6X  *  f(Xo)  +  ~|  «5X  +  ...  (3-5) 

oX  |  xo 

Using  Equation  (3—1)  (and  dropping  higher  order  terms),  this 
may  be  simplified  to 


6X  =  — I  6X  (  3—6  ) 

OX  |xo 

This  system  of  equations  is  called  the  equations  of 
variation  and  is  often  written 

<5X  =  [  A  ]  <5X  (3-7) 
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c 


■ 


where  (A),  the  variational  matrix,  is  defined  as 


[  A  ( t )  J 


XO 


(  3—8 ) 


The  elements  of  (A(t)l  may  be  found  in  Appendix  C. 

Equation  (3-7)  is  a  linear  set  of  equations  and, 
therefore,  has  a  fundamental  solution  set  of  six  independent 
solutions,  <fi  ,  whose  initial  conditions  are  4>  =  6  . 
The  Kroenecker  Delta,  6..,  equals  one  when  i  =  j  and  zero 
when  i  *  j  .  A  general  solution  to  Equation  (3—7)  may  be 
written 

6X(t)  ■  W.<t)  6X  (to)  j  =  1-6  (3-9) 

lJ  J 

In  matrix  form,  this  may  be  written 


l 


6X(t)  =  *(t,to)«5X(to)  (3-10) 

where  »(t,to),  the  state  transition  matrix,  is  a  square 
matrix  with  columns  made  of  ^  .  It  relates  the  variations 
at  time  t  to  changes  in  the  initial  conditions.  Taking  the 
derivative  of  Equation  (3-10),  substituting  into  Equation 
(3-7),  and  making  the  appropriate  cancellations  results  in 

«(t,to)  =  ( A 1# ( t , to)  (3-11) 
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Hence,  numerical  integration  of  the  equations  of  aotion 
propagates  the  reference  trajectory;  numerical  integration 
of  Equation  (3—11)  along  with  the  equations  of  motion  also 
propagates  the  displaced  trajectories  caused  by  different 
initial  conditions. 

Locating  Periodic  Orbits 

Since  information  is  now  available  on  hov  the  initial 
conditions  affect  the  orbit  trajectory,  a  method  may  be 
developed  to  determine  the  proper  initial  conditions  which 
will  produce  a  closed  orbit  for  a  given  period. 

First,  an  initial  estimate  for  the  initial  conditions 
is  needed  to  Integrate  the  equations  of  motion.  In  this 
case,  the  initial  conditions  may  be  obtained  from  circular 
orbit  calculations  (taking  into  account  the  rotating 
reference  frame.) 

Second,  the  equations  of  aotion  and  variation  will  be 
integrated  for  a  specified  amount  of  tine. 

Third,  if  the  orbit  does  not  close  after  one  period, 
then  the  Initial  conditions  should  be  altered,  and  the 
equations  reintegrated,  until  a  periodic  orbit  is  found. 
Corrections  to  the  initial  conditions  will  be  determined 
using  a  method  described  by  Wiesel  (8:117-121). 

Figure  3  shows  part  of  an  orbit  starting  on  the  y-axis 
and  continuing  for  half  a  period.  For  orbits  limited  to 
perpendicular  crossings  of  the  y  axis  in  Phobos'  orbital 
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Figure  3:  Periodic  Orbit  Symmetry 

I 


plane,  symmetry  allows  certain  conditions  to  hold  true. 
These  conditions  are  used  to  define  a  vector,  G,  such  that 

I 


G 


'  X(ti) 
Py ( t 1 ) 
P*(ti) 


O 


(3-12) 


If  G  is  evaluated  on  a  reference  trajectory  which  is  not  a 
periodic  orbit,  then  G  will  not  be  zero.  It  will  be  an 
error  vector,  e.  The  smaller  e  is,  the  closer  the 
trajectory  is  to  a  periodic  orbit.  For  orbits  that  deviate 
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fro*  the  reference  trajectory,  Xo(t),  6  may  be  rewritten 
using  a  Taylor  series  expansion  about  Xo(t) 


G|Xo(ti)  ♦  <5X(ti),  t« 


]  ■  •[' 


X(ti), 


CB16X(ti) 


+ 


where 


IB] 


OG 

3x  Xo(ti) 


’  1  0  0  0  0  ‘ 
0  0  0  1  0 
0  0  0  0  1 


If  the  deviated  orbit  is  actually  the  desired 
orbit,  then 


G  Xo(ti)  4  6X ( t* ) ,  ti  *  O 


For  the  reference  trajectory,  as  indicated. 


>|x(t»),  tij  * 


Equation  (3—13)  *ay  now  be  rewritten 


-  e  *  lBJ<5X(tt) 


(3-13) 


(3-14) 


periodic 


(3-15) 


(3-16) 


(3-17) 


What  is  needed,  however,  is  a  relationship  between  the  error 
and  the  change  in  initial  conditions.  This  may  be 


accomplished  by  using  the  state  transition  matrix.  Equation 
(3—10)  may  be  substitued  into  Equation  (3—17)  to  give 

-  e  *  (B)#(t»,to)<5X(to)  (3-18) 

Solving  for  dX ( to )  requires  the  product  (B]*(ti,to)  to  be 
Inverted;  however,  the  product  is  not  invertable  because  it 
is  a  3  x  6  matrix.  Yet,  examination  of  the  product  reveals 

'dX(to)  ' 
dY(to) 

[  B  ]$  ( ti,  to  )dX  ( to)  =  [  B  )ft  ( ti,  to)  (3-19) 

dPy(to) 

dP*(to) 

Because  of  the  way  the  problem  has  been  set  up,  the  initial 
conditions  for  X,  Py,  and  P*  are  not  allowed  to  vary,  and 
dX  =  <5Py  *  6Px  *  0.  Therefore,  the  corresponding  rows  of 
the  product  (BJ*  may  be  eliminated.  Equation  (3—19)  may  be 
rewritten 

dy (to)  ' 

l  B  )#  ( t»,  to)dX  ( to)  =  J(B]#(ti,to)  JM4  •  dZ(to)  (3-20) 

dPx(to) 
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Equation  (3—18)  may  now  be  solved  for  6X(to),  resulting  in 


'  6y (to)  ’ 

•  6Z ( to )  • 

SPk ( to ) 


£l  B  )#  ( ti,  to ) 


(3-21) 


The  initial  conditions  say  now  be  corrected  with 
Equation  (3—21),  and  the  equations  of  motion  and  variation 
reintegrated.  This  process  is  continued  until  the  elements 
of  6  are  at  an  accepted  tolerance. 


Verification  of  Equations 

Before  a  periodic  solution  may  be  found,  the  equations 
of  motion  and  variation  must  be  calculated  and  programed 
correctly. 

Equations  of  Motion  The  Hamiltonian  for  this  system  is 
not  a  function  of  time  and  should  remain  constant  through 
out  the  integration.  Therefore,  the  Hamiltonian  at  different 
points  along  the  orbit  may  be  compared.  If  they  differ, 
then  the  equations  of  motion  may  not  be  correct.  In  order 
to  locate  where  the  problem  may  lie,  different  components 
may  be  "turned  off"  (e.g.,  by  setting  the  mass  of  Mars  to 
zero)  or  they  may  be  "turned  up"  (e.g.,  by  increasing  the 
moments  of  inertia  for  Phobos).  By  varying  the  parameters 
and  observing  how  the  Hamiltonian  is  affected,  the  error  in 
the  equations  may  be  found. 
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If  the  Hamiltonian  remains  constant  during  part  of  the 
trajectory,  then  "jumps"  to  a  new  value  (and  remains 
constant  at  this  new  value),  the  trajectory  has  passed  too 
close  to  a  singularity.  The  trajectory  which  results  is  not 
reliable  and,  therefore,  the  equations  should  be 
reintegrated  with  different  initial  conditions. 

Equations  of  Variation  The  variational  matrix  may  be 
verified  by  using  its  definition 


U(t»'  *  M 


Using  an  approximation  for  the  derivative 


_  f(X  t  AX)  —  f (X) 
iaj  - 


(3-22) 


(3-23) 


where  AX  is  a  small  number.  Each  column  of  [A]  will  be 
checked  by  varying  the  corresponding  element  of  X. 

State  Transition  Matrix  The  state  transition  matrix 
may  be  verified  by  noting 


#(t,to)  = 


Mo  Xo 


(3-24) 


Using  an  approximation  for  the  derivative 


X*(t)  -  X»(t) 
t(t,to)  *  - ^ - 


(3-25) 


r 
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where  AXo  is  a  snail  change  in  the  initial  condition.  Xi(t) 
is  the  state  at  tine  t  started  fron  initial  condition  Xo, 
and  X2(t)  is  the  state  at  tine  t  started  fron  initial 
condition,  Xo  +  AXo.  Each  colunn  of  #(t,to)  will  be  checked 
by  varying  the  corresponding  elenent  of  Xo. 


( 
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IV.  Stability  and  Control 

Once  a  periodic  orbit  has  been  located,  the  next  step 
is  to  determine  if  it  is  stable.  In  other  words:  If  a 
satellite  is  placed  in  a  periodic  orbit,  will  it  remain  in 
that  orbit  (even  if  it  experiences  minor  perturbations?)  If 
the  orbit  is  inherently  unstable,  then  a  control  system 
needs  to  be  Implemented  to  ensure  stability  (2:  671-676). 

Stability  of  Orbits 

The  system  of  equations 

6X  =  ( A ( t )  1  <5X  (4-1) 

is  time  periodic  and  linear.  The  solution  for  time 
periodic,  linear  systems  was  discovered  by  Floquet.  The 
most  common  application  of  Floquet  theory  is  to  determine 
the  stability  of  periodic  orbits. 

As  discussed  previously,  a  solution  to  Equation  (4—1) 
is  given  by 


<5X(t)  =  *(t,to)6X(to) 

where  #(t,to)  is  calculated  by 


(4-2) 


#(t,to) 


lA(t) ]*(t,to) 


(4-3) 


with  the  initial  condition  for  f(t,to)  being  II],  the 
identity  matrix.  Floquet's  theorem  states  that  the  state 
transition  matrix  may  be  factored  such  that 

t(t,to)  =  lF(t)l  expUJlt)  IF(to))"1  (4-4) 

(J),  a  constant  matrix  usually  in  the  Jordan  normal  form, 
has  diagonal  entries  called  Poincar*  exponents  which  are 
labeled  u>i.  They  are  analogous  to  the  eigenvalues  of  a 
constant  coefficient  system.  The  matrix  (F(t)l  is  a  time 
periodic  matrix  with  the  same  period,  T,  as  the  original 
system.  Therefore,  IF(T)]  =  IF(0)]  and  Equation  (4-4)  may 
be  written 


#(T,0)  *  ( F(  0  )  ]  expUJJT)  I F  ( 0 ) )  ~4  (4-5) 

If  Equation  (4-5)  is  rewritten  as 

expUJJT)  =  ( F  ( 0 ) )  ~4#  ( T ,  0 )  [F(0)J  (4-6) 

it  shows  that  (F(0))  is  the  matrix  of  eigenvectors  of 
#(T,0).  The  eigenvalues,  of  t(T,0)  are  related  to  the 

Poincar*  exponents  by 
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A  =  expluf.T} 

V  V 


(4-7) 


or 


»  ( 1/T )  InA  <  4—8  ) 

t  V 

In  order  to  determine  stability.  Equation  (4-3)  is 
integrated  for  one  period  to  obtain  #(T,0).  The  Poincare 
exponents  are  then  calculated;  if  any  of  then  have  positive 
real  parts,  the  system  is  unstable. 

Calculating  t  F  ( t )  1 

The  control  system  for  this  problem  (to  be  developed 
later)  will  require  lF(t)l~4  to  be  known  at  any  point  in 
time.  Differentiating  (F(  t) )  IF(  t )  )~4  *  II)  produces 

(F(t)lCF(t))"4  +  (F(t)  JlF(t)  ]  ~4  =  10]  (4-9) 

Solving  for  IF(t)]-4 

IF(t)]"1  *  -  lF(t )  ]"4IF(t )  1  IF  ( t )  l"1  (4-10) 

An  expression  for  (F(t)]  is  now  needed.  Taking  the 
derivative  of  Equation  (4-4)  gives 

t(t,to)  *  I F ( t )  ]  etJH  IF( to)  )_1 

+  (F(t)llJl  etJn  IF(to))]'4  (4-11) 
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Substituting  Equations  (4-11)  and  (4-4)  into  Equation  (4-3) 
and  rearranging  produces 


l F( t )  )  =  lA(t))(F(t))  -  [ F ( t ) } [ 3 ]  (4-12) 

Substitution  of  Equation  (4-12)  into  (4-10)  yields  an 
expression  for  (F(t)l-1 

[  F(  t )  J_1  *  [JlIF(t)]"1  -  (F(t)  )-1lA(t)  )  (4-13) 

Integration  of  Equation  (4-12)  over  one  period  will  produce 
values  for  [F(t)l-1  for  any  time,  t.  The  initial  condition 
for  l F ( t ) )  1  is  the  inverse  of  the  matrix  of  eigenvectors  of 
•  ( T, 0  )  . 

The  Poincare  exponents  and  eigenvalues  for  the  system 
■ay  be  conplex,  but  nay  be  Manipulated  to  produce  real 
■atrices  for  both  [J1  and  (F(t)l.  (F(t)l  is  arranged  with 
the  eigenvectors,  F.,  as  its  columns.  If  F^  is  complex  (and 
hence  F  is  its  complex  conjugate),  then  the  4th  column  of 
(F(t)J  will  be  the  real  part  of  and  the  (i+l)lh  column 
will  be  the  imaginary  part  of  F^, 

[J]  will  be  arranged  in  a  block  diagonal  form.  If 
is  real,  it  will  be  a  diagonal  element  of  U).  If  ^  is 
complex  (with  being  its  complex  conjugate),  then  IJI 

will  have  a  block  diagonal  of  the  form 
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(4-14) 


I 


Re  ( u>  )  I  ■  ( u>  ) 

\  \ 

— Ia(o>  )  Re(o>  ) 

V  V  J 

Since  { F ( 0 )  ]  and  (J)  are  now  available.  Equation  (4-13) 
may  be  integrated. 

A  closed  torn  expression  is  needed  to  represent 
[F(t)]~*  so  that  its  value  at  any  tine  will  be  available. 
Since  (F(t)]-1  is  periodic,  each  of  its  elements,  f~*,  may 
be  written  as  a  Fourier  series 

f-1  =  4  C  +  C  cos  ( t )  +  C  cos  (  2t )  +  ...  +  ^  C  cos  ( nt ) 

ij  2  O  1  2  2  n 

+  S  sin(t)  +  S  sin(2t)  +  ...  +  S  sin(nt)  (4—15) 

*2  n 

where  the  coefficients  are  calculated  by 

2n-  1 

=  f~‘(la)  cos(kla),  k  =  0,1,2, ...,n 

1  =0 

(4-16) 

2  n- 1 

^  V  f74(lo()  sin(kla),  k  =  0,l,...,n-l 

*  n  /  -  '■i 

l  SO 

where  a  *  2n/2n  (1:  109).  Therefore,  numerically  finding 

2n  values  of  IF(t))-4  spaced  over  equal  time  intervals,  and 
using  those  values  in  Equations  (4-16)  and  (4-15)  will 
provide  an  expression  for  IF(t)]-1. 
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Control 


Modal  Variables  A  new  variable  is  now  introduced  by 

6X(t)  =  I F { t )  )  7><t)  (4-17) 

where  r>(t)  is  the  modal  coordinate  vector.  Taking  the 
derivative  of  the  above  expression,  and  rearranging, 
produces 

n(t)  =  l  F  ( t )  ]  £  <5X(t)  -  (F(t)ln(t)  j  (4-18) 

Substitution  of  Equation  (4-1)  and,  consequently. 
Equation  (4-17),  results  in 

r>(t )  =  (F(t)r1  J  (A(t)llF(t)!  -  [F(t)J  j  n(t)  (4-19) 
Using  Equation  (4-12)  reduces  this  to 

rj(t)  =  (J]r)(t)  (4-20) 

Therefore,  (F(t)l  is  a  periodic  transformation  which 
transforms  a  periodic  system  of  equations  to  a  constant 
coefficient  system. 

Modal  Control  Theory  The  Poincar*  exponents  for  the 
system  may  be  changed  by  adding  feedback 

6X(t)  -  (A(t)]6X(t)  +  B(t)u(t)  (4-21) 
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where  B(t)  is  a  vector  that  determines  the  distribution  of 
control.  For  this  investigation,  B(t)  is  considered 
constant  and  given  by 


B  = 


(4-22) 


which  adds  the  control  to  the  Px  direction.  Equation  (4-21) 
■ay  be  rewritten  with  modal  variables  as 

n(t)  =  U]7?(t)  +  g  ( t )  u  ( t )  (4-23) 


g(t),  the  controllability  matrix,  is 


g(t )  * l F ( t )  ] ~*B 


- i 

i « 
- 1 
14 
-  1 
•  4 
-  1 
44 
-  1 
94 
-  1 
04 


(4-24) 


where  f~*  represent  the  individual  elements  of  (F(t)]~ft. 
The  1th  mode  of  this  system  may  be  controlled  only  if  the 
ith  value  of  g(t)  is  nonzero. 


32 


For  this  system,  only  one  mode  vill  be  unstable  and 
need  to  be  controlled.  To  control  just  one  node,  u(t)  takes 


the  fora 


u ( t )  *  k(t)n(t) 


■  [°' 0 . m11'  °] 


n(  t ) 


(4-25) 


where  the  unstable  modal  variable,  r>  ( t),  is  multiplied  by 
k^t).  For  this  investigation,  k  is  considered  constant. 
Equation  (4-23)  is  nov  rewritten  as 


*  [ 


(31  +  g(t)k  T)c(t) 


(4-26) 


or  as 


nc(t)  = 


’  co 

1 

0 

0 

0  ‘ 

0 

co 

X 

0 

0 

0 

0 

kg .  +co. 

V  V 

0 

0 

0 

0 

CO 

n 

T>c(t) 


(4-27) 


where  n  (t)  denotes  the  controlled  system.  The  controlled 

C 

system  has  the  same  Poincar*  exponents  as  the  original 
system  except  for  the  ilh  exponent.  A  method  is  now  needed 
to  determine  k  in  order  to  make  co.  any  chosen  value. 

I 

The  equation  for  the  unstable  mode  is 


(4-28) 


( 
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The  use  of  an  Integrating  factor  produces  a  solution  to 
Equation  (4-28) 


r>  (t )  =  tj  (0)  exp 

cv.  cv 


[  l  h  *  ] 


(4-29) 


Since  g(t)  actually  represents  one  of  the  elements  of 
(Fit)]*1,  it  is  a  periodic  function,  and  aay  be  expressed  as 
a  Fourier  series.  It  has  a  constant  tern,  g  ,  and  a 
periodic  tern,  gip(t).  This  allovs  Equation  (4-29)  to  be 
written 


rt  At)  *  T)  (0)  exp 

cv  cv 


1  « 
[ct  +  k^ic)t]  e*P  J  kg.p(t)dt 


(4-30) 


The  second  exponential  is  periodic,  so  the  new  Poincar* 
exponent  is  the  argument  of  the  first  exponential 

«.*«.+  Kg.  ( 4—31 ) 

CV  V  7VC 


where  denotes  the  controlled  Poincare  exponent. 

Now  that  the  desired  gain  has  been  established,  the 
control  needs  to  be  iapleaented  in  the  physical  variables. 
Using  Equations  (4-24)  and  (4-17),  the  control  aay  be 
written  using  physical  coordinates  as 


(4-32) 


'c 


■ 


f 


t: 


u  ( t )  =  k  [  F  ( t ) )  ~*6X  ( t ) 

and  the  controlled  system  takes  the  £ora 

dX(t)  =  £  (A(t))  +  Bk  ( F  ( t )  1 -1  J  <5X(t)  (4-33) 

Controlling  the  State  When  the  control  system  has  been 
implemented  and  the  shift  in  the  unstable  Poincar*  exponent 
verified,  the  next  step  is  to  put  the  control  in  the 
equations  of  motion 

X(t)  =  f(X)  +  Butt)  (4-34) 


I 


Using  Equation  (4-32)  results 

X ( t )  =  f(X)  +  Bk(F(t)l" 


in 

|  X(t)  -  Xnom(t) 


} 


(4-35) 


where  Xnom(t)  is  the  trajectory  of  the  desired  periodic 
orbit.  Therefore,  if  the  current  trajectory  is  the  desired 
orbit,  there  will  be  no  control  term  added  to  the  system. 

In  summary,  controlling  an  unstable  orbit  requires  the 
following  steps: 

1.  Using  the  method  described  in  Section  III,  locate  a 
closed  orbit  for  a  specified  period,  and  determine  its 
Poincar*  exponents.  Assuming  the  orbit  is  unstable. 


i 
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the  Fourier  series  for  (Ftt)]'1  and  Xnom(t)  are  also 
needed . 

2.  Select  a  value  for  to.  and  determine  the  value  of  k 

vc 

needed  to  obtain  the  shift  in  Poincar*  exponents. 

3.  Reintegrate  the  equations  of  motion  and  variation 
using  the  controlled  system.  Equation  (4-33),  and 
verify  that  the  unstable  root  has  shifted  to  the 
desired  to.  . 

VC 

4.  The  control  is  then  Implemented  in  the  equations  of 
motion  using  Equation  (4-35)  vhere  the  effect  of  the 
control  may  be  seen  directly. 


V.  Results 


Periodic  Orbits 

The  equations  of  aotion  and  variation  were  programed 
and  verified  using  the  aethods  described  in  section  III. 
The  algorithm  for  locating  periodic  orbits  vas  then 
programed  and  debugged. 

Many  estiaates  for  initial  conditions,  both  for 
prograde  and  retrograde  orbits,  were  tried  before  the 
prograa  would  converge  to  a  set  of  initial  conditions  which 
would  produce  a  closed  orbit.  However,  the  faaily  of  orbits 
which  resulted  were  unrealizable  orbits  because  they  went 
through  Phobos. 

In  order  to  find  physically  realizable  orbits,  the 
following  steps  were  followed: 

1.  Phobos*  aoaents  of  Inertia  were  set  to  zero,  thus 
aaklng  Phobos  a  point  aass. 

2.  A  closed  orbit  which  did  not  pass  through  Phobos  was 
located . 

3.  The  aoaents  of  Inertia  were  then  increased  a  saall 
aaount . 

4.  The  algor itha  was  rerun  with  the  converged  values 
froa  Step  2  as  the  initial  conditions. 

5.  Steps  2  through  5  were  repeated  until  a  closed 
orbit  vas  located  with  the  aoaents  of  Inertial  at  their 
noraal  value. 


37 


Using  this  Method/  periodic  orbits  were  located  which  did  not 
pass  through  Phobos.  The  Initial  conditions  to  which  the 
algorithm  converged  for  three  of  these  orbits  are  shown  In 
Table  I  and  match  the  orbits  found  by  Werner  (6:23-24). 
Plgure  4  shows  Orbit  fl,  the  orbit  for  which  the  remainder  of 
the  results  are  presented. 


Table  I: 

Initial 

Conditions  for  Three 

Periodic 

Orbits 

Orbit  « 

Period 

(sec)  Y(to)  (km) 

P*(to) 

(km/sec ) 

1 

8200 

14.16398 

2.14778 

2 

8300 

14.64673 

2.14765 

3 

8400 

15.14859 

2.14754 

Stability  and  Control 

The  eigenvalues,  eigenvectors,  and  Poincar6  exponents 
were  calculated  next,  and  are  shown  in  Table  II.  The  first 
Poincar*  exponent  is  positive;  hence,  the  orbit  is  unstable. 
Therefore,  aodal  control  (described  in  Section  IV)  was  used 


to  shift  the  exponent  to  the  left  half  plane 


Calculating  t F ( t ) 3  For  this  investigation,  the 
unstable  Poincar6  exponent  is  w  ;  therefore,  only  the  first 
row  of  l F ( t )  1  ~ 1  will  need  to  be  represented  by  its  Fourier 
series.  The  Fourier  coefficients  were  calculated  using  the 
method  described  in  Section  IV,  and  are  listed  in  Appendix  D. 
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Note  that  IF(O))-1  =  IF(T)  l-1,  which  can  be  used  to  verify 
correct  prograaaing  and  integration  of  Equation  (4-13). 


Table  II:  Eigenvalues  (X),  Poincar*  Exponents  (<•>), 
and  Eigenvectors  (E)  for  Orbit  tfl 


3 . 075e  +  00  +  O.OOOe+OOi 


1 . 370e-04  +  O.OOOe+OOi 


3 . 252e-01  +  O.OOOe+OOi 


o>  *-l .  370e-04  +  O.OOOe  +  OOi 
2 


X  =  1 . 000e  +  00  +  3 . 082e-06 i 


w  *  0.000e+00  +  O.OOOe+OOi 


X4=  1 . 000e+00  -  3 . 082e-06 i 


<*>  *  0 . 000e  +  00  +  O.OOOe+OOi 


X  ■-7 . 988e-01  +  6.016e-01i 

9 


w  a  0 . 000e  +  00  +  3 . 044e-04 i 


X  *-7 . 988e-01  -  6 . 016e-01i 

9 


w  »  0 . 000e  +  00  -  3 . 044e-04 i 

9 


-1 . 296e+02 

+ 

0 . 000e+00 

i  ' 

5 . 068e  +  01 

+ 

0.000e+00 

i 

0 . 000e+00 

+ 

0.000e+00 

i 

-3 . 843e-02 

+ 

0 . 000e+00 

i 

3 . 569e-02 

+ 

0 . 000e  +  00 

i 

0 . 000e+00 

+ 

0 . 000e+00 

i 

a 

-1.314e+02 

+ 

0 . 000e+00 

i  ' 

-5 . 138e+01 

+ 

0 . 000e+00 

i 

0 . 000e+00 

+ 

0 . 000e+00 

i 

3 . 897e-02 

+ 

0 . 000e+00 

1 

3 . 619e-02 

+ 

0 . 000e  +  00 

i 

0 . 000e+00 

+ 

0 . 000e  +  00 

i 

1 . 623e-01 

— 

7 . 758e  +  03 

i  ' 

-8 . 840e-03 

- 

1 . 974e-07 

i 

0 . 000e+00 

+ 

0 . 000e+00 

i 

2.434e-06 

+ 

0 . 000e+00 

i 

-1 . 231e-04 

+ 

5 . 884e+00 

i 

0 . 000e+00 

+ 

0 . 000e  +  00 

i 

a 

1 . 623e-01 

+ 

7 . 758e  +  03 

i  ‘ 

-8 . 840e-03 

+ 

1 . 97  4e-07 

i 

0 . 000e+00 

♦ 

0 . 000e+00 

i 

2 . 434e-06 

♦ 

0 . 000e  +  00 

i 

-1 . 231e-04 

- 

5 . 884e+00 

i 

0 . 000e+00 

+ 

0 . 000e  +  00 

i 

0 . 000e+00 

+ 

0 . 000e+00 

i 

0 . 000e+00 

+ 

0 . 000e  +  00 

1 

-3 . 559e+01 

+ 

2 . 726e  +  02 

i 

0 . 000e+00 

+ 

0 . 000e  +  00 

1 

0 . 000e+00 

+ 

0 . 000e+00 

1 

3 . 339e-01 

+ 

4 . 360e-02 

i 

0 . 000e+00 

+ 

0 . 000e  +  00 

i 

0 . 000e+00 

+ 

0 . 000e+00 

i 

-3 . 559e+01 

- 

2 . 726e+02 

i 

0 . 000e+00 

+ 

0 . 000e  +  00 

i 

0 . 000e+00 

+ 

0 . 000e  +  00 

i 

3 . 339e-01 

• 

“ 

4 . 360e-02 

i 

Modal  Control  Examination  of  Equation  (4—31)  shows 


w  =  a>.  +  kg  ( 5—3 ) 

a  v  vc 

where,  in  this  case,  the  unstable  exponent  Is 

=  1 . 370e-04  (5-4) 

Shifting  the  unstable  exponent  to  a  desired  value  of 

«ci  =  -2 . 000e-04  (5-5) 

results  In  the  equation 

-2 . 000e-04  =  1 . 370e-04  +  kg  (5-6) 

1C 

g  is  the  constant  part  (the  C©  Fourier  coefficient)  of  f_1 

1C  14 

(see  Equation  (4—24))  and  is 

g  *  -14.3244  (5-7) 

ic 

Equation  (5—6)  nay  now  be  solved  for  k,  resulting  in 

k  =  2.3526e-05  (5-8) 

I 
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Equation  (4—33)  any  now  be  integrated  to  verify  the  shift  in 
the  Poincare  exponent.  The  Poincar*  exponents  for  the 
controlled  systea  are  listed  in  Table  III/  and  aatch 
predicted  results. 

Table  III:  Poincar*  Exponents  for  the  Controlled  Systea 
to  =  -2 . 000e-04  +  0.000e  +  00  i 

t 

to  *  -1 . 370e-04  +  0 . 000e  +  00  i 
2 

to  =  0.000e  +  00  +  0 . 000e  +  00  i 

a 

to  =  0 . 000e  +  00  +  0 . 000e  +  00  i 

4 

to  *  0 . 000e+00  ♦  3 . 044e-04  i 

9 

to  *  0 . 000e  +  00  -  3 . 044e-04  i 

<j 

I  - 

State  Control  Equation  (4—35)  aay  now  be  used  to  add 
the  control  to  the  state.  Various  initial  conditions  were 
used  to  test  how  the  control  system  reacted  to  changes  in 
the  noainal  trajectory.  Figure  5  shows  the  amount  of 
control  thrust  needed  to  maintain  the  nominal  orbit  when  the 
initial  conditions  for  X  were  displaced  from  the  nominal 
value  by  .5,  10,  and  50  meters.  The  periodic  oscillations 
which  appear  on  the  graph  are  caused  by  the  truncation  of  the 
Fourier  series  for  (F(t)l  1  and  X(t).  The  sign  change  for 
the  initial  control  value  aay  also  be  explained  by  the 
Fourier  series  truncations.  These  phenomenon  aay  be  more 
readily  observed  by  plotting  the  control  output  for  a  nominal 
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trajectory,  shown  in  Figure  6.  The  small  error  in  the 
Fourier  series  representations  of  [F(t))~*  and  X(t)  produce  a 
control  output,  which  starts  as  a  positive  value  and  then 
converges  to  the  same  periodic  oscillations  as  do  the  curves 
in  Figure  5. 

When  the  initial  conditions  are  displaced  beyond  the 
linear  control  region,  the  control  becomes  unstable.  This 
is  shown  in  Figure  7.  For  the  X  direction,  the  instability 
occurs  when  the  displacement  exceeds  5  kilometers. 

The  control  thrust  was  also  plotted  for  initial 
displacements  in  Y ,  Px,  and  Py  for  both  the  stable  and 
unstable  regions.  These  graphs  are  similar  to  Figures  5  and 
7,  and  are  shown  in  Appendix  F.  Table  IV  lists  the 
displacements  in  X,  Y,  Px,  Py  which  produce  unstable  control 
outputs  when  the  control  thrust  is  in  the  Px  direction. 

Table  IV:  Displacements,  Which  Produce  Unstable 
Control  in  Px  Direction 

AX  :  5  km  APx  :  . 5  cm/s ec 

AY  :  10  m  APy  :  1  m/sec 


It  should  also  be  noted  that  displacements  in  Z  or  Pz 
are  not  controllable  because  f~*(t)  and  f~*(t)  are  zero  for 

II  10 

all  time  (See  Equation  (4—35)).  However,  the  orbit  is 
already  stable  in  these  two  modes,  as  investigation  of  the 
Poincar6  exponents  ,  <*>#  and  o0,  will  reveal. 
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The  control  output  vas  also  investigated  vhen  the 
thrust  vas  applied  in  the  Py  direction.  This  is 

accomplished  by  setting  ( B 1  to 


The  gain  also  needs  to  be  recalculated  using  Equation  (5—3) 

where  g  is  now  the  constant  part  of  £_1.  The  new  gain  is 
’ic  is  ’ 


k  =  4 . 0188e-Q5 


(5-10) 


The  resulting  control  plots  were  similar  to  the  plots 
presented  above.  However,  the  region  of  control  instability 
changed,  as  shown  in  Table  VI. 


Table  v:  Displacements.  Which  Produce  Unstable 
Control  in  Py  Direction 

AX  :  10  m  APx  :  .5  cm/sec 

AY  :  10  a  APy  :  .5  cm/sec 
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VI.  Conclusions  and  Recommendations 


Periodic  orbits  about  Phobos*  equatorial  plane  have 
been  located  and  discovered  to  be  inherently  unstable. 
These  orbits  may  be  stabilized  by  adding  a  control  system 
which  will  shift  the  unstable  Poincar*  to  the  left  hand 
plane.  When  the  control  system  is  implemented  in  the  Pm 
direction,  the  amount  of  displacement  away  from  the  nominal 
trajectory  which  can  be  controlled  successfully  is  greater 
than  if  the  control  were  in  the  Py  direction.  Therefore, 
putting  the  control  system  in  the  P*  direction  would  be  more 
beneficial . 

A  possible  extension  of  this  research  would  be  to 
investigate  other  types  of  orbits  around  Phobos  to  determine 
their  stability  and,  if  needed,  controllability. 


Appendix  A:  Problea  Paraneters 

The  following  values  were  used  when  Integrating  the 
equations  of  notion  and  variation  (6:35). 

Phobos  Axis  Lengths 
x  =  10.6  ka 
y  *  13.5  ka 
z  =  9 . 4  ka 

(These  values  are  used  to  coapute  the  aoaents  of 
Inertia . ) 

Phobos'  Orbital  Radius 
D  =  9378  km 

Phobos  Rotation  Rate 

O  =  2.28  x  10'4  sec"* 

GMrhob  »  6.6  x  10  4  kmVsec* 

GMmot.  =  42828.32  kmVsec* 
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Appendix  B:  Equations  of  Motion 

The  equations  of  motion  were  derived  in  Section  II,  and 
listed  here  for  reference. 


X  =  2S-  *  Px  +  f>( y  -  D) 
<7Px 


(2 


y 


9H 

apx 


Py  -  OX 


(2 


2 


OH 

OPx 


P* 


(2- 


Px 


—  *  OPy 
9X  y 

GMmor* 

X  _ 

3XY 

■ 

3R*X 

GMphobX 

D 

D* 

D* 

2D4 

R- 

+  _  25*  { 5x*  -  2R*MI  -  2Ixx) 

4R9  4R?  _ 


+  5Y* ( I  -  Iyy)  +  5Z*(I 


I«) 


1 


are 

36a) 

36b) 

36c) 


(2-36d) 


A  _  GMmor. 

Py  *  “  5y  “  ~CPx - D — 


GHpKobY  3GIY 


.[  1  _  4Y_  3(YR*~  DR 
"I  D  D*  2D* 

l.m\  5X*(  I  -  21. 
S  4R  [ 


*-  2DY*] 


4  ( 5Y*—  2R*)(I  -  Iyy )  +  5Z*(I  -  Ik) 


( 2— 36e ) 


P*  =  - 


GM  mar  i 

5T~ 


r.  Z  _  3YZ  _  3R2Z  1  _  GMpKobZ 

D*  D*  2D*  J  R* 

—  -  —  5X*(  I  -  2Ixx ) 
t  4R7 


4  5Y*(  I  -  Iyy)  4  ( 5Z*-  2R*)(I  -  Ik) 


( 2— 36f ) 


Appendix  C:  Equations  of  Variation 


The  equations  of  variation  are  defined  as  (See  Section  III) 


1Alt"  -  3*,. 


Therefore,  (A(t)J  has  the  fora 


(3-8) 


( A  ( t )  ) 


■  0 

O 

0 

1 

0 

o  ' 

-o 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

l 

A«t 

A«2 

A4I 

0 

n 

0 

Asi 

AS2 

Asa 

-n 

0 

0 

Am 

AS2 

Asa 

0 

0 

0 

(C-l) 


The  eleaents  for  the  lover  left  block  are  listed  below. 


GMmots  3Y  .  3(2X*  +  R*) 

A«t  *  J.  T 

D  D  2D 


GMrhob ( R2—  3X2)  3GI(R*~  5X*) 

RS  4R7 


♦  ^2R4—  25X*R*t  35X4J  £l  -  2I*x} 

+  5^7X2-  R*J£y*(I  -  2Iyy )  +  Z*(I  -  21**) J 


(C— 2 ) 
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Asz  = 


Ass  = 


GMmots 

D* 


4  +  3  (  2Y*  +  R1  -  6DY  ) 

2D* 


GMf hob ( R2—  3YZ)  +  3GI (R2-  5Y2) 


4R 


3G 


4R 


^2R4—  25Y*R*+  35Y4j  -  2Iy 


(7Y2-  R*)(> 


+  5|7Y‘-  R‘||X*(I  -  2Ikx)  +  Z*(I  - 


\ 


3GMxarsZ ( Y  -  D)  .  3GMshobYZ 
-  +  - 


15GYZI 

4R7 


+ 


15GYZ 

4RP 


2Iyy] 


♦  ^7Z*—  2R*J  £l  -  2I**J  +  7X*£l  - 


Am  =  Ass 


Appendix  D:  Fourier  Coefficients  for  [ F ( t ) 1 


-t 


Table  D.l:  Fourier  Coefficients  for  f*1 

it 


n 

Cn 

Sn 

0 

-3.874357399e-03 

•  •  • 

1 

-8. 78375710 4e -05 

-1.328999018e-03 

2 

-1. 24367240 le -03 

-3. 07078402 6e -03 

3 

-5. 80432828 5e -04 

-1.177822744e-03 

4 

-1. 86947339 8e -04 

-3. 49588202 3e -04 

5 

-7. 48229624 4e -05 

-1. 309451163c -04 

6 

-2.281366727e-05 

-3. 14144494 8e -05 

7 

3.760515418e-06 

1. 26686884 3e -05 

8 

7. 385247274c -06 

1. 42140935 le -05 

9 

2. 96886255 5e -06 

4 . 5891 2 720 2e -06 

10 

7.251091917e-08 

-3.653094817e-07 

11 

-4. 41126039 2e -07 

-8. 74432404 5e -07 

12 

-1. 88786515 le -07 

-2. 95582183 8e -07 

Table  D.2:  Fourier  Coefficients  for  f 

12 


n 

Cn 

Sn 

0 

2 . 117141101e-03 

•  •  • 

1 

-2.690985867e-03 

2 . 112497840e-04 

2 

-3.351871379e-03 

1 . 321959626e-03 

3 

-1 . 21 8 5 326 50e -03 

5 . 963226336e-04 

4 

-3. 53171811 3e -04 

1. 88611209 3e -04 

5 

-1. 28258994 4e -04 

7 . 345863939e-05 

6 

-2. 97010412 8e -05 

2. 18513730 5e -05 

7 

1 . 293546998e-05 

-3. 95612694 8e -06 

8 

1 . 405638000e-05 

-7. 31622172 4e -06 

9 

4 . 491670339e-06 

-2. 91376560 5e -06 

10 

-3. 78763 500 8e -07 

-6. 09 0322784c -08 

11 

-8.666678737e-07 

4.377670807e-07 

12 

-2. 91509860 3e -07 

1.865378556e-07 

13 

1. 02051464 9e -08 

9. 93711341 le -09 

14 

4. 37125612 Oe -08 

-2. 134858245c -08 

Table  D.3:  Fourier  Coefficients  for  f-1 


n 

Cn 

Sn 

0 

0. 00000000 Oe +00 

•  »  • 

1 

0. 000000000e+00 

0. 00000000 Oe +  00 

2 

0. 00000000 Oe +00 

0. 00000000 Oe +  00 

3 

0. 00000000 Oe +00 

0.000000000e+00 

4 

0.000000000e+00 

0. 00000000 Oe +00 

5 

0. 00000000 Oe +00 

0. 00000000 Oe +  00 

6 

0. 00000000 Oe +00 

0.000000000e+00 

7 

0.000000000e+00 

0.000000000e+00 

8 

0.00 000000 Oe +00 

0. 00000000 Oe +  00 

9 

0.000000000e+00 

0.000000000e+00 

10 

0. 00000000 Oe +00 

0. 00000000 Oe +00 

11 

0.000000000e+00 

0. 00000000 Oe +00 

12 

0. 00000000 Oe +00 

0. 00000000 Oe +00 

13 

0.000000000e+00 

0 . 000000000e+00 

14 

0. 00000000 Oe +00 

0. 00000000 Oe +00 

Table  D.4:  Fourier  Coefficients  for  f 

14 


n 

Cn 

Sn 

0 

1. 43244116 2e +  01 

•  •  • 

1 

2. 89299384 Oe +  00 

-6. 00043048 le -01 

2 

2 . 455985030e+00 

7 . 010133811e-01 

3 

-5.866175316e-01 

2 . 420906168e-01 

4 

-1.262025226e-01 

5 . 984262457e-02 

5 

-3.703386356e-02 

1 . 933474339e-02 

6 

-7. 32877157 5e -03 

4 . 979233267e-03 

7 

2.485946774e-03 

-6. 67375013 3e -04 

8 

2. 431689468e-03 

-1. 19386891 5e -03 

9 

6. 9906 36 18 le -04 

-4. 27949682 le -04 

10 

-5. 30337985 le -05 

-1.397491958e-05 

11 

-1. 06792534 3e -04 

5. 20813563 8e -05 

12 

-3. 48374733 8e -05 

2. 32263263 le -05 

13 

8. 46551147 3e -07 

-2.656660668e-06 

14 

3.572232102e-06 

1.003728450e-06 

59 


Table  D.5:  Fourier  Coefficients  for  f-1 

ts 


n  Cn  Sn 


0 

-8.385572830e+00 

•  •  • 

1 

-6.634642137e-01 

4. 25412768  Oe  4-00 

2 

7.328663541e-01 

2 . 618125520e+00 

3 

2.474867601e-01 

6. 03024073 7e -01 

4 

6 . 029949451e-02 

1 . 273106157e-01 

5 

1.9024088856-02 

3 . 635979353e-02 

6 

4 . 792466561e-03 

6 . 966502661e-03 

7 

-7.016951413e-04 

-2. 53681294 5e -03 

8 

-1.1856 50 449e-03 

-2.412165999e-03 

9 

-4.255307518e-04 

-6. 81983429 8e -04 

10 

-9. 52549177 8e -06 

5.256663100e-05 

11 

5.285072527e-05 

1. 07123976 3e -04 

12 

2.124386919e-05 

3. 21 58 4526 5e -05 

13 

1.579965902e-06 

-2.073022030e-06 

14 

-2. 8 56 46752 6e -06 

-6. 32985304 6e -06 
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Table 

D.6:  Fourier  Coefficients  for  f-t 

id 

*'  ’l 

4 

a 

n 

Cn 

Sn 

0 

0. 00000000 0e +00 

•  •  • 

pfc 

1 

0. 00000000 Oe +00 

0.000000000e+00 

i 

2 

0. 00000000 Oe +  00 

0. 00000000 Oe +  00 

3 

0. 00000000 Oe +00 

0 . 000000000e  +  00 

4 

0.000000000e+00 

0. 00000000 Oe +  00 

5 

0. 00000000 Oe +00 

0.000000000e+00 

« 

6 

0. 00000000 Oe +00 

0. 00000000 Oe +  00 

7 

0.000000000e+00 

0 . 000000000e  +  00 

8 

0. 00000000 Oe +  00 

0. 00000000 Oe +  00 

9 

0.000000000e+00 

0.000000000e+00 

4 

10 

0 . 000000000e  +  00 

0, 00000000 Oe +  00 

* 

11 

0. 00000000 Oe +  00 

0. 00000000 Oe +00 

12 

0. 00000000 Oe +  00 

0.000000000e+00 

■ 

13 

0.000000000e+00 

0. 00000000 Oe +00 

i 

14 

0. 00000000 Oe +  00 

0.000000000e+00 

ft 
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Appendix  E:  Fourier  Coefficients  for  X(t) 


Table  E.l:  Fourier  Coefficients  for  X 


n 

Cn 

Sn 

0 

.  88  3488e -11 

•  •  • 

1 

.  140728e-09 

.  152364e+02 

2 

.  135547e-10 

.  6  49805e  +  00 

3 

.  100078e-ll 

. 544625e-01 

4 

. 871430e-12 

. 949703e-02 

5 

- .  127151e-12 

.14921 3e -02 

6 

. 598870e-13 

- . 765536e-03 

7 

- . 247681e-13 

- . 561069e-03 

8 

- .  545792e-13 

- . 127291e-03 

9 

.  769867e-14 

. 166649e-04 

10 

- . 256899e-13 

. 192629e-04 

11 

- . 505583e-13 

. 467848e-05 

12 

- .  24  3318e-13 

- . 436199e-06 

13 

- . 347421e-14 

- . 611304e-06 

14 

-  . 630388e-13 

- . 151512e-06 
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Table  E.2:  Fourier  Coefficients  for  Y 
n  Cn  Sn 


0 

.  479096e+01 

•  •  • 

1 

. 182127e+02 

- . 169381e-09 

2 

.  677786e*00 

-.118613e-10 

3 

. 551567e-01 

.  212204e-ll 

4 

.  928735e-02 

.  481067e-ll 

5 

.  139793e-02 

.  670131e-ll 

6 

.  774109e-03 

.  821775e-ll 

7 

, 554801e-03 

. 965711e-ll 

8 

.  124648e-03 

.110613e~10 

9 

.  168713e~04 

. 124592e-10 

10 

.  190833e-04 

.  138822e-10 

11 

. 460648e-05 

.  152965e-10 

12 

.  441411e-06 

.  167062e-10 

13 

.  606409e-06 

. 181132e-10 

14 

.  149646e-06 

.  I95201e-10 

Table  E.3:  Fourier  Coefficients  for  Z 


i 


n 


Cn 


Sn 


0 

.  OOOOOOefOO 

«  •  • 

1 

. OOOOOOefOO 

.  000000e  +  00 

2 

. 000000e+00 

.OOOOOOefOO 

3 

. 000000e+00 

.  000000e  +  00 

4 

. OOOOOOefOO 

.  OOOOOOefOO 

5 

. 000000e+00 

.OOOOOOefOO 

6 

.OOOOOOefOO 

.  OOOOOOefOO 

7 

. 000000e+00 

.  000000e  +  00 

8 

. 000000e+00 

.OOOOOOefOO 

9 

•OOOOOOefOO 

.  000000e  +  00 

10 

.OOOOOOefOO 

.OOOOOOefOO 

11 

.OOOOOOefOO 

.OOOOOOefOO 

12 

. OOOOOOefOO 

.OOOOOOefOO 

13 

. OOOOOOefOO 

.OOOOOOefOO 

14 

. OOOOOOefOO 

.OOOOOOefOO 

Table  E.4:  Fourier  Coefficients  for  Px 


n 

Cr» 

Sn 

0 

.  213928e+01 

•  •  • 

1 

.  752224e-02 

. 14274 5e -12 

2 

.  841281e-03 

. 407597e-12 

3 

. 112619e-03 

. 635699e-12 

4 

.  269906e-04 

. 847254e-12 

5 

.  539791e-05 

. 106258e-ll 

6 

- .  334303e-05 

. 127447e-ll 

7 

-.288292e-05 

. 148699e-ll 

e 

- . 751878e-06 

. 170020e-ll 

9 

.  111064e-06 

. 191237e-ll 

10 

.  143236e-06 

. 212447e-ll 

11 

. 383699e-07 

. 23373 5e -11 

12 

- .  392334e-08 

. 254974e -11 

13 

- « 596416e-08 

. 27S121e-ll 

14 

- .  1604  32e-08 

. 297555e-ll 

Table  E.S:  Fourier  Coefficients  for  Py 


n 


Cn 


Sn 


0 

. 249991e-14 

•  •  • 

1 

- . 978226e-13 

- . 104815e-01 

2 

- . 183950e-13 

- . 890540e -03 

3 

- . 363208e-14 

- . 114373e-03 

4 

- . 125753e-14 

- . 263001e-04 

5 

- . 253473e-15 

- . 501555e-05 

6 

. 117735e-15 

.338438e-05 

7 

. 148175e-15 

. 284786e-05 

8 

. 247384e-16 

. 735059e-06 

9 

- . 599974e-16 

- . 112549e-06 

10 

- . 651474e-16 

- . 141833e-06 

11 

- . 537480e-16 

37760 5e -07 

12 

- . 552417e-16 

. 395848e-08 

13 

- . 593738e-16 

. 590029e-08 

14 

- . 492293e-16 

. 156984e-08 
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Table  E.6:  Fourier  Coefficients  for  Pz 
n  Cn  Sn 


0 

. 000000e+00 

•  •  • 

1 

. 000000e+00 

. 000000e+00 

2 

. 000000e+00 

.  000000e  +  00 

3 

. 000000e+00 

. 000000e+00 

4 

. 000000e+00 

.  000000e  +  00 

5 

. 000000e+00 

.  000000e  +  00 

6 

. 000000e+00 

.  000000e  +  00 

7 

. 000000e+00 

. 000000e+00 

8 

. 000000e+00 

.  OOOOOOe  +  00 

9 

. 000000e+00 

. 000000e+00 

10 

. 000000e+00 

.  000000e  +  00 

11 

. 000000e+00 

. 000000e+00 

12 

. 000000e+00 

. 000000e+00 

13 

.000000e+00 

. 000000e+00 

14 

. 000000e+00 

.  000000e  +  00 

20.00 


000 1 


000 1 


o 

Fq 

t  CO 


p-TT 


(zoas/uj>!  0(0  Lx)  |OJ}uoo 
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Figure  F.4:  Nonlinear  Control  for  Displacements  in  Px 
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